Numerical Simulation of Three-Dimensional Dendritic Growth 
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Dendritic crystal growth in a pure undercooled melt is simulated quantitatively in three dimensions 
using a phase-field approach. The full non-axisymmetric morphology of the steady-state dendrite 
tip and a* are determined as a function of anisotropy for a crystal with a cubic symmetry. Results 
are compared to experiment and used to critically test solvability theory. (June, 1996). 
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Dendritic growth has been a central problem in pattern 
formation |l| and metallurgy for many years. Consid- 
erable theoretical progress has been achieved which has 
led to the development of solvability theory to determine 
the steady-state operating state of the dendrite tip (i.e. 
tip velocity V and tip radius p) [|l]j|~|^]. 

Paradoxically, our understanding of pattern selection 
in three dimensions (3D) has remained uncertain, espe- 
cially with regards to experiment. This is due to the fact 
that it has so far remained too difficult to simulate re- 
liably the equations of dendritic growth in 3D even on 
current supercomputers. This, in turn, has prevented: 
(i) to test whether the global attractor of the growth dy- 
namics is indeed the steady-state needle crystal predicted 
by solvability theory, and (ii) to test the predictions of 
this theory which are themselves only approximate in 3D 
|^-^| . Consequently, it has remained unclear whether ex- 
isting disagreements between solvability theory and ex- 
periment p-|Tl|l are due to the inapplicability of this the- 
ory, to the approximate nature of its predictions in 3D, 
or to some missing physics in the starting equations. 

Here we report the results of 3D simulations of den- 
dritic growth in a pure undercooled melt. These simu- 
lations are made possible by using a recently developed 
phase-field approach which renders fully quantitative 3D 
computations accessible for the first time jl2]. This ap- 
proach also makes it possible to model the experimentally 
relevant low velocity limit where the solid- liquid interface 
can be assumed to relax instantaneously to local ther- 
modynamic equilibrium (i.e. where kinetic effects at the 
interface are negligibly small) . Our numerical results are 
applied to critically test the applicability of solvability 
theory in 3D, to make comparison with experiment, and 
to characterize the full 3D dendrite tip morphology. 

The equations of the 3D free-boundary problem are 
given by 
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where T is the temperature field, Tj is the interface tem- 
perature, Tm is the melting temperature, D is the diffu- 



sivity, L is the latent heat of melting, c p is the specific 
heat, v n is the normal velocity of the interface, and 



7( n ) = 7o (1 ~ 3e 4 ) 
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is the surface energy for a crystal chosen to have a cubic 
symmetry. In addition, 9i are the angles between the 
normal n and the two local principal directions on the 
boundary and Ri are the principal radii of curvature. 
Growth is controlled by the undercooling A = (Tm — 
Too) / (L / c p ) where T^ is the initial melt temperature. 

We use three key ingredients to solve the above equa- 
tions. Firstly, we avoid the usual difficulties of tracking a 
sharp boundary by using a phase-field approach @ Q . 
In this approach, the two-phase system is described by a 
phcnomcnological free-energy 



T= / dr[M/ 2 (n)|V^| 2 + /(V,w)] 
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where u = (T — Tm)/(L/c p ) is the dimensionless tem- 
perature field. We have used here the form f(tp, u) = 
-VJ 2 /2 + V> 4 /4 + Xuip(l- 2vj 2 /2, + V> 4 /5) with minima 
at -0 = — 1 an d ip = +1 that correspond to the liq- 
uid and solid phases, respectively. The anisotropic sur- 
face energy defined by Eq. ^ is recovered by choosing 
W(n) = W 7(n)/7o, with n = Vip/\Vip\, which is the 
direct 3D generalization of the way anisotropy has been 
previously included in 2D ]T^-p^|]. The phase-field, ip, 
varies between its two minima values across an interfa- 
cial region of width Wo, thereby distinguishing between 
phases without front-tracking. Its dynamics, defined by 
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is then coupled to that of the diffusion field 
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in such a way that the equations for the two fields re- 
duce to those of the original solidification equations in 
the sharp-interface limit where Wq is small compared to 
the principal radii of curvature of the boundary. 

Secondly, and most importantly, we exploit the re- 
sults of a recent analysis [O to relate the phase-field 
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FIG. 1. Results of 3D phase-field simulation on a 
300 x 300 x 300 cubic lattice for e 4 = 0.047 which shows 
dendrite tips growing along the principal < 100 > directions. 
The simulations was performed in the first octant (a;, y, z > 0) 
with a spherical nucleus centered at the origin as initial con- 
dition. The solid-liquid boundary shown here corresponds 
to the tp = surface reconstructed by reflection about the 
x — y — z — planes. The structure is seen from an angle 
where all six < 100 > directions are visible. 

model parameters to those of the sharp-interface equa- 
tions. As demonstrated in 2D |jl2), this analysis makes 
it possible: (i) to perform computations with a smaller 
capillary length, do = jqTmCp/ L , and (ii) to choose A 
and the function r(n) in such a way that we obtain a 
Gibbs-Thomson condition without interface kinetics in 
the sharp-interface limit Jljj . That is we recover exactly 
Eq. without an additional velocity-dependent term in 
this limit. Since the computation time scales as (do /Wo) 5 
in 3D, this first property dramatically enhances the com- 
putational efficiency of the phase-field approach and al- 
lows us to model dendritic growth quantitatively, as op- 
posed to non-quantitatively in 3D as previously Jl4| . Fur- 
thermore, it permits us to model an undercooling range 
where the dendrite tip Peclet number P = pV/2D is 
small enough to compare a* = 2Ddo/p 2 V to its mea- 
sured small undercooling values ]9|,|Io|]. 

Lastly, we are able to resolve numerically very small 
anisotropics by incorporating quantitatively the contri- 
bution of the grid anisotropy. For a given value of £4 used 
in W(n), we compute the equilibrium shape produced by 
the phase-field model. The grid-corrected anisotropy is 
then defined as that value of £4 in 7(11) which produces 
a 3D equilibrium shape that matches exactly that of the 
phase-field model. A procedure that will be described 
elsewhere was also developed to obtain a grid-corrected 
r(n). Numerical tests were performed in 2D to check 
that this procedure yields accurate values of a* over the 
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FIG. 2. Plot of a* vs £ 4 for A = 0.45 showing the results 
of phase-field simulations compared to the approximate pre- 
dictions of the numerical and linear solvability theories. Also 
plotted are the small A experimental values of a* for SCN [9] 
and PVA [10] using the £4 measurements of Ref. [11]. 

entire range of £4 investigated here in 3D. 

Eqs. were simulated explicitly with a grid spacing 
ranging between 0.6 and 0.8 with Wo = 1, D ranging 
between 0.5 and 4, a time step ranging between 0.01 
and 0.08, and a constant undercooling A = 0.45. The 
anisotropy was varied from £4 = 0.0066 to £4 = 0.047 
which spans most of the range of experimental interest. 
Fig. 1 shows a typical 3D dendrite morphology resulting 
from the growth of a small spherical seed. Quantita- 
tive results which pertain to the steady-state operating 
state and morphology of the dendrite tip are presented 
in Table 1 and Figs. 2-3. These were obtained with long 
simulations that focused on the steady-state growth of a 
single dendrite tip along the z-axis. Several runs were 
performed to check that the values in Table 1 are inde- 
pendent of do /Wo within an accuracy of 10%. 

In order to test solvability theory, we have computed 
independently the values of P and a* predicted by the 
numerical solution of the steady-state growth equations 
using the standard axisymmetric approximation where 
the surface energy and steady-state shape are assumed 
to be independent of the polar angle <f> in the x-y plane 
perpendicular to the growth axis [^ . This is the same cal- 
culation performed by Kessler and Levine in Ref. j(| and 
we have checked that our boundary integral code repro- 
duces their results within one percent. For completeness, 
we also report the predictions of the linear solvability 
theory of Barbieri and Langer Q] which uses a shape lin- 
earized around the Ivantsov paraboloid of revolution p"§| ] 
and the same axisymmetric approximation. 

As shown in Table 1, the numerical solvability the- 
ory yields a* values which are systematically lower than 
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TABLE I. Result of phase-field simulations on a 
200 x 200 x 400 cubic lattice compared to the results of the nu- 
merical and linear solvability theories for A = 0.45. A and a 
characterize the amplitude of the four-fold symmetry compo- 
nent of the tip morphology in simulations. Typical runs took 
60-140 CPU hours on a DEC-ALPHA 3000-700 workstation 
and shorter times on a CRAY-YMP and a CRAY-T3D. 



C4 



0.0066 
0.0123 
0.0171 
0.0265 
0.0369 
0.0470 



Phase-Field Simulations 



.4 



0.426 0.015 

0.360 0.036 

0.312 0.063 

0.236 0.16 

0.159 0.47 

0.093 1.72 



1.78 
1.73 
1.68 
1.62 
1.57 
1.54 



13.0 
11.3 
10.1 
7.8 
5.7 
4.0 



Solvability theory 



Numerical 



0.418 0.0128 

0.367 0.0329 

0.324 0.0578 

0.247 0.142 

0.172 0.365 

0.109 1.037 



Linear 



Plv 



0.471 
0.471 
0.471 
0.471 
0.471 
0.471 



0.0116 
0.0250 
0.0367 
0.0602 
0.0890 
0.1240 



the phase-field results, but still reasonably close for small 
anisotropy. The predicted values of a* , however, start to 
become significantly inaccurate for anisotropics greater 
than about 3%, although predictions of the Peclet num- 
ber remain relatively accurate. Most likely, this does not 
indicate a breakdown of solvability theory, but of the ax- 
isymmetric approximation. This conclusion is supported 
by the fact that the four-fold deviation from a shape of 
revolution increases in magnitude with anisotropy (as de- 
scribed below) and, hence, can affect the selection. It is 
also supported by the fact that we do not observe any 
sidebranching |l|,^,^9| without adding noise to the phase- 
field equations. Hence our simulations rule out the possi- 
bility of a dynamical attractor other than a steady-state 
needle crystal. The linear theory is seen to breakdown 
at much smaller anisotropy. This is because it assumes 
that the steady-state shape remains close to the Ivantsov 
paraboloid of revolution. Table 1 shows that the actual 
Peclet number already starts deviating significantly from 
its Ivantsov value, Pi v ^M, at small anisotropy. 

The steady-state morphology of the dendrite tip was 
analyzed using the Fourier decomposition 



b, z) — A n (z) cos 4n0 
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where r is the radial distance from the z-axis. Both r 
and z are measured in units of p with the tip at z = 0. 
This decomposition has the advantage that it is gen- 
eral and does not presupposes a particular analytical 
form to fit the tip shape. The function Aq(z) is the ax- 
isymmetric contribution to this shape. It approaches a 
paraboloid of revolution, Aq(z) — ~2z, for small \z\ and 
departs from this shape with increasing \z\ as shown in 
Fig. 3, this departure being more pronounced at larger 
anisotropy as one would expect. More important is the 
non-axisymmctric departure from this shape of revolu- 
tion which is contained in the higher modes, A n (z) for 
n > 1. The amplitude of the first four-fold symmetry 



mode turns out to be much larger than all the other 
modes (Fig. 3) and to be well-described by the power law 
A\(z) = A \z\ a , which appears as a remarkably straight 
line on a log-log plot of Ai(z) vs \z\ up to a distance of a 
few p behind the tip. Furthermore, the values of A and a 
in Table 1 clearly show that the amplitude of the four- fold 
symmetry mode is sensitively dependent on anisotropy 
which is a qualitatively novel aspect of our results. In 
contrast, the linear solvability calculation of Ben Amar 
and Brener Q predicts that for small £4 the tip morphol- 
ogy should be independent of anisotropy, with A^ 1 = 11 
and a = 2 for small \z\. The values of A and a listed 
in Table 1 indicate that this prediction is most likely 
only valid for values of £4 which lie outside the range of 
experimental interest. This is also consistent with the 
fact that their calculation is only valid in the limit where 
a* ~ e^ 4 , and that this 7/4 power law scaling is not yet 
attained at the smallest computed anisotropy in Fig. 2. 

For £4 = 0.0066, which lies inside the range of uncer- 
tainty of the measured value £4 = 0.0055 ± 0.0015 for 
SCN our simulations yield a value of a* w 0.015 
with P — 0.426. This value is only 20 % smaller than 
the value a* = 0.0192 measured by Huang and Glicks- 
man (9) for this material. Most of this discrepancy can 
be accounted for by the finite Peclet number correction 
which can be estimated, using the numerical solvability 
code, to decrease a* by about 15% from its zero Peclet 
number limiting value. Therefore we obtain a reasonably 
good agreement with experiment for SCN within the ex- 
isting uncertainty in the measured value of anisotropy. 
Pivalic acid |l^] (PVA), however, remains problematic as 
seen in Fig. 2 and a closer examination of kinetic ef- 
fects for this material could potentially resolve this large 
discrepancy. Finally, Maurer et al. |^0| have found that 
the tip morphology of NH4Br dendrites is indeed well- 
described by a single cos 4</> mode, while LaCombe et al. 
have found that for SCN dendrites more modes seem nec- 
essary to fit the tip morphology |^lj . The origin of this 
difference, which is potentially due to noise amplification, 
also remains to be understood. 

In conclusion, we have demonstrated that quantita- 
tive modeling of 3D dendritic growth is possible using a 
recently developed phase-field methodology [|l2| . Our re- 
sults are important in that they leave little doubt about 
the conceptual validity of solvability theory in 3D and 
the importance of crystalline anisotropy, and yield results 
which are consistent with experiment, at least for SCN 
where kinetic effects are believed to be small. At present, 
the applicability of solvability theory — in terms of mak- 
ing accurate predictions — remains mainly limited by the 
axisymmctric approximation. Our results have shown 
that the three-dimensional morphology of the dendrite 
tip is dominated by a four-fold component whose ampli- 
tude depends sensitively on anisotropy. This prediction 
should be testable experimentally. A host of other mi- 
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crostructural pattern formation issues, such as the de- 
termination of the steady-state shape further behind the 
tip and sidebranching, are now amenable to quantitative 
study using the present computational approach. 
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FIG. 3. Steady-state tip morphology for C4 = 0.0123 
shown in: (a) the 4> = 0° (solid line) and <j> = 45° (dash-dotted 
line) planes, and (b) (100) planes equally spaced along z by 
one tip radius. Length is measured in units of p. The Ivantsov 
parabola z — — r /2 (dash-line) and a circle of unit radius are 
superimposed in (a) . The Fourier amplitudes at two tip radii 
behind the tip in (b) are: Ai(-2) = 0.29, A 2 (-2) = 0.020, 
and As(— 2) = 0.0022, illustrating that the tip morphology is 
dominated by the four-fold symmetry mode. 
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